home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / liboctave / dRowVector.cc < prev    next >
C/C++ Source or Header  |  1997-03-07  |  7KB  |  390 lines
  1. // RowVector manipulations.
  2. /*
  3.  
  4. Copyright (C) 1996 John W. Eaton
  5.  
  6. This file is part of Octave.
  7.  
  8. Octave is free software; you can redistribute it and/or modify it
  9. under the terms of the GNU General Public License as published by the
  10. Free Software Foundation; either version 2, or (at your option) any
  11. later version.
  12.  
  13. Octave is distributed in the hope that it will be useful, but WITHOUT
  14. ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  15. FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  16. for more details.
  17.  
  18. You should have received a copy of the GNU General Public License
  19. along with Octave; see the file COPYING.  If not, write to the Free
  20. Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
  21.  
  22. */
  23.  
  24. #if defined (__GNUG__)
  25. #pragma implementation
  26. #endif
  27.  
  28. #ifdef HAVE_CONFIG_H
  29. #include <config.h>
  30. #endif
  31.  
  32. #include <iostream.h>
  33.  
  34. #include "f77-fcn.h"
  35. #include "lo-error.h"
  36. #include "mx-base.h"
  37. #include "mx-inlines.cc"
  38. #include "oct-cmplx.h"
  39.  
  40. // Fortran functions we call.
  41.  
  42. extern "C"
  43. {
  44.   int F77_FCN (dgemv, DGEMV) (const char*, const int&, const int&,
  45.                   const double&, const double*,
  46.                   const int&, const double*, const int&,
  47.                   const double&, double*, const int&,
  48.                   long);
  49.  
  50.   double F77_FCN (ddot, DDOT) (const int&, const double*, const int&,
  51.                    const double*, const int&);
  52. }
  53.  
  54. // Row Vector class.
  55.  
  56. bool
  57. RowVector::operator == (const RowVector& a) const
  58. {
  59.   int len = length ();
  60.   if (len != a.length ())
  61.     return 0;
  62.   return equal (data (), a.data (), len);
  63. }
  64.  
  65. bool
  66. RowVector::operator != (const RowVector& a) const
  67. {
  68.   return !(*this == a);
  69. }
  70.  
  71. RowVector&
  72. RowVector::insert (const RowVector& a, int c)
  73. {
  74.   int a_len = a.length ();
  75.   if (c < 0 || c + a_len > length ())
  76.     {
  77.       (*current_liboctave_error_handler) ("range error for insert");
  78.       return *this;
  79.     }
  80.  
  81.   for (int i = 0; i < a_len; i++)
  82.     elem (c+i) = a.elem (i);
  83.  
  84.   return *this;
  85. }
  86.  
  87. RowVector&
  88. RowVector::fill (double val)
  89. {
  90.   int len = length ();
  91.   if (len > 0)
  92.     for (int i = 0; i < len; i++)
  93.       elem (i) = val;
  94.   return *this;
  95. }
  96.  
  97. RowVector&
  98. RowVector::fill (double val, int c1, int c2)
  99. {
  100.   int len = length ();
  101.   if (c1 < 0 || c2 < 0 || c1 >= len || c2 >= len)
  102.     {
  103.       (*current_liboctave_error_handler) ("range error for fill");
  104.       return *this;
  105.     }
  106.  
  107.   if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
  108.  
  109.   for (int i = c1; i <= c2; i++)
  110.     elem (i) = val;
  111.  
  112.   return *this;
  113. }
  114.  
  115. RowVector
  116. RowVector::append (const RowVector& a) const
  117. {
  118.   int len = length ();
  119.   int nc_insert = len;
  120.   RowVector retval (len + a.length ());
  121.   retval.insert (*this, 0);
  122.   retval.insert (a, nc_insert);
  123.   return retval;
  124. }
  125.  
  126. ColumnVector
  127. RowVector::transpose (void) const
  128. {
  129.   return ColumnVector (*this);
  130. }
  131.  
  132. RowVector
  133. real (const ComplexRowVector& a)
  134. {
  135.   int a_len = a.length ();
  136.   RowVector retval;
  137.   if (a_len > 0)
  138.     retval = RowVector (real_dup (a.data (), a_len), a_len);
  139.   return retval;
  140. }
  141.  
  142. RowVector
  143. imag (const ComplexRowVector& a)
  144. {
  145.   int a_len = a.length ();
  146.   RowVector retval;
  147.   if (a_len > 0)
  148.     retval = RowVector (imag_dup (a.data (), a_len), a_len);
  149.   return retval;
  150. }
  151.  
  152. RowVector
  153. RowVector::extract (int c1, int c2) const
  154. {
  155.   if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
  156.  
  157.   int new_c = c2 - c1 + 1;
  158.  
  159.   RowVector result (new_c);
  160.  
  161.   for (int i = 0; i < new_c; i++)
  162.     result.elem (i) = elem (c1+i);
  163.  
  164.   return result;
  165. }
  166.  
  167. // row vector by row vector -> row vector operations
  168.  
  169. RowVector&
  170. RowVector::operator += (const RowVector& a)
  171. {
  172.   int len = length ();
  173.  
  174.   int a_len = a.length ();
  175.  
  176.   if (len != a_len)
  177.     {
  178.       gripe_nonconformant ("operator +=", len, a_len);
  179.       return *this;
  180.     }
  181.  
  182.   if (len == 0)
  183.     return *this;
  184.  
  185.   double *d = fortran_vec (); // Ensures only one reference to my privates!
  186.  
  187.   add2 (d, a.data (), len);
  188.   return *this;
  189. }
  190.  
  191. RowVector&
  192. RowVector::operator -= (const RowVector& a)
  193. {
  194.   int len = length ();
  195.  
  196.   int a_len = a.length ();
  197.  
  198.   if (len != a_len)
  199.     {
  200.       gripe_nonconformant ("operator -=", len, a_len);
  201.       return *this;
  202.     }
  203.  
  204.   if (len == 0)
  205.     return *this;
  206.  
  207.   double *d = fortran_vec (); // Ensures only one reference to my privates!
  208.  
  209.   subtract2 (d, a.data (), len);
  210.   return *this;
  211. }
  212.  
  213. // row vector by matrix -> row vector
  214.  
  215. RowVector
  216. operator * (const RowVector& v, const Matrix& a)
  217. {
  218.   RowVector retval;
  219.  
  220.   int len = v.length ();
  221.  
  222.   int a_nr = a.rows ();
  223.   int a_nc = a.cols ();
  224.  
  225.   if (a_nr != len)
  226.     gripe_nonconformant ("operator *", 1, len, a_nr, a_nc);
  227.   else
  228.     {
  229.       int a_nr = a.rows ();
  230.       int a_nc = a.cols ();
  231.  
  232.       if (len == 0)
  233.     retval.resize (a_nc, 0.0);
  234.       else
  235.     {
  236.       // Transpose A to form A'*x == (x'*A)'
  237.  
  238.       int ld = a_nr;
  239.  
  240.       retval.resize (a_nc);
  241.       double *y = retval.fortran_vec ();
  242.  
  243.       F77_XFCN (dgemv, DGEMV, ("T", a_nr, a_nc, 1.0, a.data (),
  244.                    ld, v.data (), 1, 0.0, y, 1, 1L));
  245.  
  246.       if (f77_exception_encountered)
  247.         (*current_liboctave_error_handler)
  248.           ("unrecoverable error in dgemv");
  249.     }
  250.     }
  251.  
  252.   return retval;
  253. }
  254.  
  255. // other operations
  256.  
  257. RowVector
  258. RowVector::map (d_d_Mapper f) const
  259. {
  260.   RowVector b (*this);
  261.   return b.apply (f);
  262. }
  263.  
  264. RowVector&
  265. RowVector::apply (d_d_Mapper f)
  266. {
  267.   double *d = fortran_vec (); // Ensures only one reference to my privates!
  268.  
  269.   for (int i = 0; i < length (); i++)
  270.     d[i] = f (d[i]);
  271.  
  272.   return *this;
  273. }
  274.  
  275. double
  276. RowVector::min (void) const
  277. {
  278.   int len = length ();
  279.   if (len == 0)
  280.     return 0;
  281.  
  282.   double res = elem (0);
  283.  
  284.   for (int i = 1; i < len; i++)
  285.     if (elem (i) < res)
  286.       res = elem (i);
  287.  
  288.   return res;
  289. }
  290.  
  291. double
  292. RowVector::max (void) const
  293. {
  294.   int len = length ();
  295.   if (len == 0)
  296.     return 0;
  297.  
  298.   double res = elem (0);
  299.  
  300.   for (int i = 1; i < len; i++)
  301.     if (elem (i) > res)
  302.       res = elem (i);
  303.  
  304.   return res;
  305. }
  306.  
  307. ostream&
  308. operator << (ostream& os, const RowVector& a)
  309. {
  310. //  int field_width = os.precision () + 7;
  311.  
  312.   for (int i = 0; i < a.length (); i++)
  313.     os << " " /* setw (field_width) */ << a.elem (i);
  314.   return os;
  315. }
  316.  
  317. istream&
  318. operator >> (istream& is, RowVector& a)
  319. {
  320.   int len = a.length();
  321.  
  322.   if (len < 1)
  323.     is.clear (ios::badbit);
  324.   else
  325.     {
  326.       double tmp;
  327.       for (int i = 0; i < len; i++)
  328.         {
  329.           is >> tmp;
  330.           if (is)
  331.             a.elem (i) = tmp;
  332.           else
  333.             break;
  334.         }
  335.     }
  336.   return is;
  337. }
  338.  
  339. // other operations
  340.  
  341. RowVector
  342. linspace (double x1, double x2, int n)
  343. {
  344.   RowVector retval;
  345.  
  346.   if (n > 0)
  347.     {
  348.       retval.resize (n);
  349.       double delta = (x2 - x1) / (n - 1);
  350.       retval.elem (0) = x1;
  351.       for (int i = 1; i < n-1; i++)
  352.     retval.elem (i) = x1 + i * delta;
  353.       retval.elem (n-1) = x2;
  354.     }
  355.  
  356.   return retval;
  357. }
  358.  
  359. // row vector by column vector -> scalar
  360.  
  361. double
  362. operator * (const RowVector& v, const ColumnVector& a)
  363. {
  364.   double retval = 0.0;
  365.  
  366.   int len = v.length ();
  367.  
  368.   int a_len = a.length ();
  369.  
  370.   if (len != a_len)
  371.     gripe_nonconformant ("operator *", len, a_len);
  372.   else if (len != 0)
  373.     retval = F77_FCN (ddot, DDOT) (len, v.data (), 1, a.data (), 1);
  374.  
  375.   return retval;
  376. }
  377.  
  378. Complex
  379. operator * (const RowVector& v, const ComplexColumnVector& a)
  380. {
  381.   ComplexRowVector tmp (v);
  382.   return tmp * a;
  383. }
  384.  
  385. /*
  386. ;;; Local Variables: ***
  387. ;;; mode: C++ ***
  388. ;;; End: ***
  389. */
  390.